# Replication Code for Effort vs. Accuracy
# Author: Marika Landau-Wells (mlw@berkeley.edu)


# Table of Contents #
# 1. Import & Tidy Data
# 2. PCA
# 3. Ground-Truth Hypervolumes (x 4)
# 4. Random Responses Distributions - DO NOT RE-RUN
# 5. Mentalizer Performance
# 6. Between-Condition Tests
# 7. Logistic Regressions 
# 8. Tests vs. Random Response Benchmarks
# 9. Main Text Figures: 2-5 
# 10. Main Text Tables: 1, 2
# 11. Appendix Tables: 1, 3, 4

# Clear the environment
gc()
rm(list = ls())

# Originally created with R version 3.6.2 (2019-12-12)
# Note that changes to R's randomization procedures mean that 
# running this code using later versions of R will not generate 
# exact replications of the hypervolumes and associated analyses

# For reproducibility, run R 3.6.2 and use groundhog to load packages
library(groundhog)

pkgs = c("dplyr", "stringr", "hypervolume", "ggplot2", "cocor", "stargazer")
groundhog.library(pkgs, "2019-12-12")

# To load independently, use these versions:
#library(dplyr) #v 1.0.2
#library(ggplot2) #v 3.2.1
#library(stringr) #v 1.4.0
#library(hypervolume) #v 2.0.12
#library(cocor) #v1.1-3
#library(stargazer) #v 5.2.2

# Set the working directory
#setwd()

#################################################
# 1. Import & Tidy Data ####

## A. Main Experiment Dataset ####

# N = 839
d.expt <- readRDS(file = "MainExpt_Clean_FINAL.rds") 

# Save variable groups
demo.vars <- c("ResponseID", "Gender", "Age", "SelfID", "Edu",
               "PID_7", "Ideo_7")
cond.vars <- c("Condition", "Type", "Item", "PT_order", "Danger")
thr.vars <- c("PhysHarm_self", "PhysHarm_loved", "Loss_econ", "Loss_rights", 
              "Loss_status", "Infection", "Purity_self", "Purity_group",
              "Environ")
emo.vars <- c("Fear", "Anger", "Disgust", "Happy", "Sad", "Pity", "Resentment", 
              "Contempt","Pride", "Anxiety")

## B. Subset by Condition (8 Subsets) ####

### Issue: Immigration ####
d.imm <- filter(d.expt, Item == "Illegalimm") #408
# Self/Other Partitions
d.imm.self <- filter(d.imm, Type == "Self") #197
d.imm.oth <- filter(d.imm, Type == "Other") #211

# Block Order Partitions
# N corresponds to Figure 1A
d.imm.s.t1 <- filter(d.imm.self, PT_order == 1) #99
d.imm.s.e1 <- filter(d.imm.self, PT_order == 2) #98
d.imm.o.t1 <- filter(d.imm.oth, PT_order == 1) #108
d.imm.o.e1 <- filter(d.imm.oth, PT_order == 2) #103

### Issue: Climate Change ####
d.cch <- filter(d.expt, Item == "Climate") #431
# Self/Other Partitions
d.cch.self <- filter(d.cch, Type == "Self") #217
d.cch.oth <- filter(d.cch, Type == "Other") #214

# Block Order Partitions
# N corresponds to Figure 1A
d.cch.s.t1 <- filter(d.cch.self, PT_order == 1) #111
d.cch.s.e1 <- filter(d.cch.self, PT_order == 2) #106
d.cch.o.t1 <- filter(d.cch.oth, PT_order == 1) #108
d.cch.o.e1 <- filter(d.cch.oth, PT_order == 2) #106

## C. Define Ground-Truth #### 
d.imm.s.t1ref <- filter(d.imm.s.t1, Danger >= 50) #57
d.imm.s.e1ref <- filter(d.imm.s.e1, Danger >= 50) #62

d.cch.s.t1ref <- filter(d.cch.s.t1, Danger >= 50) #79
d.cch.s.e1ref <- filter(d.cch.s.e1, Danger >= 50) #89

#################################################

#################################################
# 2. PCA ####
## A. Inputs for PCA ####
### Ground Truth ####
# Threat Perception 
imm.t1.ref <- select(d.imm.s.t1ref, one_of(thr.vars)) #57
cch.t1.ref <- select(d.cch.s.t1ref, one_of(thr.vars)) #79

# Emotion Response 
imm.e1.ref <- select(d.imm.s.e1ref, one_of(emo.vars)) #62
cch.e1.ref <- select(d.cch.s.e1ref, one_of(emo.vars)) #89

### Mentalizers ####
# Threat-First: Threat Rating
imm.t1.o <- select(d.imm.o.t1, ResponseID, Danger, one_of(thr.vars)) #108
cch.t1.o <- select(d.cch.o.t1, ResponseID, Danger, one_of(thr.vars)) #108

# Threat-First: Emotion Rating
imm.e2.o <- select(d.imm.o.t1, ResponseID, Danger, one_of(emo.vars)) #108
cch.e2.o <- select(d.cch.o.t1, ResponseID, Danger, one_of(emo.vars)) #108

# Emotion-First: Threat Rating
imm.t2.o <- select(d.imm.o.e1, ResponseID, Danger, one_of(thr.vars)) #103
cch.t2.o <- select(d.cch.o.e1, ResponseID, Danger, one_of(thr.vars)) #106

# Emotion-First: Emotion Rating
imm.e1.o <- select(d.imm.o.e1, ResponseID, Danger, one_of(emo.vars)) #103
cch.e1.o <- select(d.cch.o.e1, ResponseID, Danger, one_of(emo.vars)) #106


# Performed mostly to deal with correlations in the data structure
# 3 PCs don't account for equal amounts of variance for imm and cch
# but always > 75% 

## B. Ground Truth PCAs: II ####
# Threat Perception (Threats-First Condition):
imm.t1.mat <- data.matrix(imm.t1.ref)
imm.t1.pca <- prcomp(imm.t1.mat, center = T, scale = F)
imm.t1.vexp = imm.t1.pca$sdev^2 / sum(imm.t1.pca$sdev^2)
sum(imm.t1.vexp[1:3]) 
# 3 PCs: 0.8515294
# dim(imm.t1.pca$x) 
# 57 x 9

# Emotion Response (Emotions-First Condition):
imm.e1.mat <- data.matrix(imm.e1.ref)
imm.e1.pca <- prcomp(imm.e1.mat, center = T, scale = F)
imm.e1.vexp = imm.e1.pca$sdev^2 / sum(imm.e1.pca$sdev^2)
sum(imm.e1.vexp[1:3]) 
# 3 PCs: 0.7585869
# dim(imm.e1.pca$x) 
# 62 x 10

## C. Ground Truth PCAs: CC ####
# Threat Perception (Threats-First Condition):
cch.t1.mat <- data.matrix(cch.t1.ref)
cch.t1.pca <- prcomp(cch.t1.mat, center = T, scale = F)
cch.t1.vexp = cch.t1.pca$sdev^2 / sum(cch.t1.pca$sdev^2)
sum(cch.t1.vexp[1:3]) 
# 3 PCs: 0.7864986
# dim(cch.t1.pca$x) 
# 79  9

# Emotion Response (Emotions-First Condition):
cch.e1.mat <- data.matrix(cch.e1.ref)
cch.e1.pca <- prcomp(cch.e1.mat, center = T, scale = F)
cch.e1.vexp = cch.e1.pca$sdev^2 / sum(cch.e1.pca$sdev^2)
sum(cch.e1.vexp[1:3]) 
# 3 PCs: 0.7739873
# dim(cch.e1.pca$x) 
# 89  10

## D. Mentalizer PCAs: II ####

### Threat Estimation ####
# Threats-First Condition
imm.o.t1.pca <- predict(imm.t1.pca, data.matrix(imm.t1.o[,-c(1:2)]))
#dim(imm.o.t1.pca)
#108 x 9

# Emotion-First Condition
imm.o.t2.pca <- predict(imm.t1.pca, data.matrix(imm.t2.o[,-c(1:2)]))
#dim(imm.o.t2.pca)
#103 x 9

### Emotion Understanding ####
# Threats-First Condition
imm.o.e2.pca <- predict(imm.e1.pca, data.matrix(imm.e2.o[,-c(1:2)]))
#dim(imm.o.e2.pca)
#108 x 10

# Emotion-First Condition
imm.o.e1.pca <- predict(imm.e1.pca, data.matrix(imm.e1.o[,-c(1:2)]))
#dim(imm.o.e1.pca)
#103 x 10

## E. Mentalizer PCAs: CC ####
### Threat Estimation ####
# Threats-First Condition
cch.o.t1.pca <- predict(cch.t1.pca, data.matrix(cch.t1.o[,-c(1:2)]))
#dim(cch.o.t1.pca)
#108 x 9

# Emotion-First Condition
cch.o.t2.pca <- predict(cch.t1.pca, data.matrix(cch.t2.o[,-c(1:2)]))
#dim(cch.o.t2.pca)
#106 x 9

### Emotion Understanding ####
# Threats-First Condition
cch.o.e2.pca <- predict(cch.e1.pca, data.matrix(cch.e2.o[,-c(1:2)]))
#dim(cch.o.e2.pca)
#108 x 10

# Emotions-First Condition
cch.o.e1.pca <- predict(cch.e1.pca, data.matrix(cch.e1.o[,-c(1:2)]))
#dim(cch.o.e1.pca)
#106 x 10

#################################################

#################################################
# 3. Ground-Truth Hypervolumes ####

## A. Inputs for HVs ####
# Specify the number of PCs on which to base the HV
pc.use = 3

# Pick the Ground-Truth data:
imm.t1.hv.df <- data.frame(imm.t1.pca$x[, 1:pc.use])
imm.e1.hv.df <- data.frame(imm.e1.pca$x[, 1:pc.use])
cch.t1.hv.df <- data.frame(cch.t1.pca$x[, 1:pc.use])
cch.e1.hv.df <- data.frame(cch.e1.pca$x[, 1:pc.use])


## B. Build HVs ####
# HV construction relies on random sampling
set.seed(1234)

# Use SVM with defaults
# Ignore error about SD because 
# preserving that property is intentional and scales are already consistent
imm.t1.hv <- hypervolume(imm.t1.hv.df, 
                         name = "imm_t1",
                         method="svm",
                         svm.nu = 0.01, 
                         svm.gamma = 0.5,
                         scale.factor = 1,
                         verbose = F)

imm.e1.hv <- hypervolume(imm.e1.hv.df, 
                         name = "imm_e1",
                         method="svm",
                         svm.nu = 0.01, 
                         svm.gamma = 0.5,
                         scale.factor = 1,
                         verbose = F)

cch.t1.hv <- hypervolume(cch.t1.hv.df, 
                         name = "cch_t1",
                         method="svm",
                         svm.nu = 0.01, 
                         svm.gamma = 0.5,
                         scale.factor = 1,
                         verbose = F)

cch.e1.hv <- hypervolume(cch.e1.hv.df, 
                         name = "cch_e1",
                         method="svm",
                         svm.nu = 0.01, 
                         svm.gamma = 0.5,
                         scale.factor = 1,
                         verbose = F)

## C. HV Stats ####
# If volumes are different from those listed, stats won't replicate
imm.t1.hv@Volume #1519265
imm.e1.hv@Volume #2117455
cch.t1.hv@Volume #1580007
cch.e1.hv@Volume #2120507

# If Number of Random Points differs, stats won't replicate
dim(imm.t1.hv@RandomPoints)
#18940     3
dim(imm.e1.hv@RandomPoints)
#21775     3
dim(cch.t1.hv@RandomPoints)
#15732      3
dim(cch.e1.hv@RandomPoints)
#14568     3
#################################################

#################################################
# 4. Random Response Distributions - DO NOT RUN ####
# Note: This code documents how simulated data were generated
# However, the HV function for distances doesn't replicate perfectly
# even using set.seed()
# For exact replication, use the provided files

# ## A. Helper Functions ####
# ### Function for Generating Random Ratings ####
# Generates a list of random draw Other-Ratings
# with same dimensions as specified via x; n datasets is length of list
rand.oth.fun <- function(x, n){
  # Note: seed set outside this function to avoid same results for both simulations
  # x is the oth.mat or df, just for dims
  # n is the number of samples to draw
  # all possible ratings
  ratings.vec <- c(0:100)
  # count of rows and columns
  ratings.n <- dim(x)[1]
  features.n <- dim(x)[2]
  # results list
  out.ls <- list()
  # create n datasets
  for(i in 1:n){
    out.mat <- t(replicate(ratings.n,
                           sample(ratings.vec,
                                  features.n,
                                  replace = T),
                           simplify = T))
    colnames(out.mat) <- colnames(x)
    out.ls[[i]] <- out.mat
  }
  return(out.ls)
}

# NOTE: the random draw does not need to differ
# between issues since the space of points is the same (101^9 or 101^10)
# But sample sizes are slightly different, so each Issue gets its own
# Randomly drawn set

# ### Function for Collecting Results from Random Ratings #### 
# (1) Inclusion Test
# (2) Euclidean distances
# (3) Accuracy + distances
hypv.rand.fun2 <- function(x, v){
  # x = a list of pca-transformed guess matrices, pc-use applied
  # v = original hypervolume
  # Note: pcs to use is value defined above
  # lists to hold distance vectors
  guess.ls <- list()
  for (i in 1:length(x)){ #n samples
    # pull matrix
    use.mat <- x[[i]]
    # inclusion test for each point as vector (T/F)
    inctest.vec <- hypervolume_inclusion_test(v,
                                              points=use.mat,
                                              reduction.factor=1,
                                              fast.or.accurate = "accurate",
                                              verbose = F)
    # iterate through rows to get distances
    dist.vec <- c()
    set.seed(1234)
    for (j in 1:length(use.mat[,1])){
      # for distance, make each row a hypervol (list)
      orow <- matrix(use.mat[j,], ncol = pc.use, nrow=1, byrow=T)
      oexb <- expectation_ball(input = orow, num.samples = 5)
      # get the distance to the original
      # set.seed(1234)
      dist.vec[j] <- hypervolume_distance(hv1 = v,
                                          hv2 = oexb,
                                          type = "minimum",
                                          num.points.max = 500,
                                          check.memory = F)
      
      rm(orow, oexb)
    }
    # Using inclusion test results and dist.vec:
    # Invert the inclusion test (0 = inside)
    inctest.miss <- as.numeric(!inctest.vec)
    # take the mean of the distances (break apart if desired)
    guess.vec <- inctest.miss * dist.vec
    guess.ls[[i]] <- data.frame(Obs = 1:length(use.mat[,1]),
                                Inc_test = as.numeric(inctest.vec),
                                Euc_dist = dist.vec,
                                Guess_dist = guess.vec)
  }
  return(guess.ls)
}



# ## B. Generate Random Responses ####
#
# ### Draw Datasets (N=500) ####
# # Note: if set.seed(1234) is run within the loop, this produces
# # nearly identical results, wrapped differently due to number of dims
# # setting seed outside function provides althernative draw
# 
# set.seed(1234)
# t1.imm500.rand <- rand.oth.fun(x=imm.t1.o, n=500)
# t1.cch500.rand <- rand.oth.fun(x=cch.t1.o, n=500)
# #head(t1.imm500.rand[[1]])
# #head(t1.cch500.rand[[1]])
# 
# set.seed(5678)
# e1.imm500.rand <- rand.oth.fun(x=imm.e1.o, n=500)
# e1.cch500.rand <- rand.oth.fun(x=cch.e1.o, n=500)
# #head(e1.imm500.rand[[1]])
# #head(e1.cch500.rand[[1]])
# 
# ## C. PCA on Random Responses ####
# ### Threat Estimation ####
# t1.imm.rand.pca <- lapply(t1.imm500.rand,
#                       function(x)
#                         predict(imm.t1.pca, x))
# t1.cch.rand.pca <- lapply(t1.cch500.rand,
#                           function(x)
#                             predict(cch.t1.pca, x))
# 
# t1.imm.rand.pcs <- lapply(t1.imm.rand.pca,
#                       function(x)
#                         matrix(x[,1:pc.use],
#                                ncol=pc.use,
#                                byrow=T))
# t1.cch.rand.pcs <- lapply(t1.cch.rand.pca,
#                           function(x)
#                             matrix(x[,1:pc.use],
#                                    ncol=pc.use,
#                                    byrow=T))
# 
# ### Emotion Understanding ####
# e1.imm.rand.pca <- lapply(e1.imm500.rand,
#                       function(x)
#                         predict(imm.e1.pca, x))
# e1.cch.rand.pca <- lapply(e1.cch500.rand,
#                           function(x)
#                             predict(cch.e1.pca, x))
# e1.imm.rand.pcs <- lapply(e1.imm.rand.pca,
#                       function(x)
#                         matrix(x[,1:pc.use],
#                                ncol=pc.use,
#                                byrow=T))
# e1.cch.rand.pcs <- lapply(e1.cch.rand.pca,
#                           function(x)
#                             matrix(x[,1:pc.use],
#                                    ncol=pc.use,
#                                    byrow=T))
# 
# ## D. Measure Binary Accuracy and Miss Distances for Random Responses ####
# ### Threat Estimation ####
# # Immigration
# system.time(
#   imm.t1.rand.dist <- hypv.rand.fun2(x = t1.imm.rand.pcs,
#                                      v = imm.t1.hv)
# )
# # Approx. 3 minutes
# #imm.t1.rand.dist[[2]]
# 
# # Climate
# system.time(
#   cch.t1.rand.dist <- hypv.rand.fun2(x = t1.cch.rand.pcs,
#                                      v = cch.t1.hv)
# )
# # Approx. 2 minutes
# #cch.t1.rand.dist[[2]]
# 
# ### Emotion Understanding ####
# # Immigration
# system.time(
#   imm.e1.rand.dist <- hypv.rand.fun2(x = e1.imm.rand.pcs,
#                                      v = imm.e1.hv)
# )
# # Approx. 3 min
# 
# # Climate
# system.time(
#   cch.e1.rand.dist <- hypv.rand.fun2(x = e1.cch.rand.pcs,
#                                      v = cch.e1.hv)
# )
# # Approx. 2 min
# 
# # Save all of the random distance objects
# rand.hv.dist <- list(imm.t1.rand.dist, cch.t1.rand.dist,
#                      imm.e1.rand.dist, cch.e1.rand.dist)
# #save(rand.hv.dist, file = "AccDist_500Random_Draws_240719.RDS")
# 
#################################################

#################################################
# 5. Mentalizer Performance ####

## A. Binary Accuracy ####
### Threat Estimation, Threats-First Conditions ####
imm.o.t1.inc <- hypervolume_inclusion_test(imm.t1.hv, 
                                           points=imm.o.t1.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(imm.o.t1.inc)/length(imm.o.t1.inc) 
#0.712963
# 108

cch.o.t1.inc <- hypervolume_inclusion_test(cch.t1.hv, 
                                           points=cch.o.t1.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(cch.o.t1.inc)/length(cch.o.t1.inc) 
#0.7962963

### Threat Estimation, Emotions-First Conditions ####
imm.o.t2.inc <- hypervolume_inclusion_test(imm.t1.hv, 
                                           points=imm.o.t2.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(imm.o.t2.inc)/length(imm.o.t2.inc) 
#0.776699
# 103

cch.o.t2.inc <- hypervolume_inclusion_test(cch.t1.hv, 
                                           points=cch.o.t2.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(cch.o.t2.inc)/length(cch.o.t2.inc) 
#0.8113208


### Emotion Understanding, Emotions-First Conditions ####
imm.o.e1.inc <- hypervolume_inclusion_test(imm.e1.hv, 
                                           points=imm.o.e1.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(imm.o.e1.inc)/length(imm.o.e1.inc) 
#0.6504854

cch.o.e1.inc <- hypervolume_inclusion_test(cch.e1.hv, 
                                           points=cch.o.e1.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(cch.o.e1.inc)/length(cch.o.e1.inc) 
# 0.7641509

### Emotion Understanding, Threats-First Conditions ####
imm.o.e2.inc <- hypervolume_inclusion_test(imm.e1.hv, 
                                           points=imm.o.e2.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(imm.o.e2.inc)/length(imm.o.e2.inc) 
#0.6388889

cch.o.e2.inc <- hypervolume_inclusion_test(cch.e1.hv, 
                                           points=cch.o.e2.pca[,1:pc.use], 
                                           reduction.factor=1,
                                           fast.or.accurate = "accurate",
                                           verbose = F)
sum(cch.o.e2.inc)/length(cch.o.e2.inc) 
#0.8703704


## B. Miss Distances ####

# NOTE: Calculations of distances reported in the publication did not 
# use set.seed() consistently.  This script produces the paper's 
# results exactly.  If set.seed() is uncommented in all eight distance
# calculations, distances (and associated statistical tests) will differ slightly.

### Issue: II ####
# Threat Estimation, Threats-First Condition
imm.o.t1.d <- c()
for (i in 1:length(imm.o.t1.pca[,1])){ 
  # pull row for sub
  o.guess <- imm.o.t1.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  set.seed(1234)
  imm.o.t1.d[i] <- hypervolume_distance(hv1 = imm.t1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
imm.t1.miss <- !imm.o.t1.inc #inverts the inclusion test
imm.t1.missdist <- imm.t1.miss * imm.o.t1.d

# Threat Estimation, Emotions-First Condition
imm.o.t2.d <- c()
for (i in 1:length(imm.o.t2.pca[,1])){ 
  # pull row for sub
  o.guess <- imm.o.t2.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  #set.seed(1234)
  imm.o.t2.d[i] <- hypervolume_distance(hv1 = imm.t1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
imm.t2.miss <- !imm.o.t2.inc #inverts the inclusion test
imm.t2.missdist <- imm.t2.miss * imm.o.t2.d

# Emotion Understanding, Threats-First Condition
imm.o.e2.d <- c()
for (i in 1:length(imm.o.e2.pca[,1])){ 
  # pull row for sub
  o.guess <- imm.o.e2.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  set.seed(1234)
  imm.o.e2.d[i] <- hypervolume_distance(hv1 = imm.e1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
imm.e2.miss <- !imm.o.e2.inc #inverts the inclusion test
imm.e2.missdist <- imm.e2.miss * imm.o.e2.d

# Emotion Understanding, Emotions-First Condition
imm.o.e1.d <- c()
for (i in 1:length(imm.o.e1.pca[,1])){ 
  # pull row for sub
  o.guess <- imm.o.e1.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  #set.seed(1234)
  imm.o.e1.d[i] <- hypervolume_distance(hv1 = imm.e1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
imm.e1.miss <- !imm.o.e1.inc #inverts the inclusion test
imm.e1.missdist <- imm.e1.miss * imm.o.e1.d



### Issue: CC ####
# Threat Estimation, Threats-First Condition
cch.o.t1.d <- c()
for (i in 1:length(cch.o.t1.pca[,1])){ 
  # pull row for sub
  o.guess <- cch.o.t1.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  #set.seed(1234)
  cch.o.t1.d[i] <- hypervolume_distance(hv1 = cch.t1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

cch.t1.miss <- !cch.o.t1.inc #inverts the inclusion test
cch.t1.missdist <- cch.t1.miss * cch.o.t1.d

# Threat Estimation, Emotions-First Condition
cch.o.t2.d <- c()
for (i in 1:length(cch.o.t2.pca[,1])){ 
  # pull row for sub
  o.guess <- cch.o.t2.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  #set.seed(1234)
  cch.o.t2.d[i] <- hypervolume_distance(hv1 = cch.t1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
cch.t2.miss <- !cch.o.t2.inc #inverts the inclusion test
cch.t2.missdist <- cch.t2.miss * cch.o.t2.d

# Emotion Understanding, Emotions-First Condition
cch.o.e1.d <- c()
for (i in 1:length(cch.o.e1.pca[,1])){ 
  # pull row for sub
  o.guess <- cch.o.e1.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  set.seed(1234)
  cch.o.e1.d[i] <- hypervolume_distance(hv1 = cch.e1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
cch.e1.miss <- !cch.o.e1.inc #inverts the inclusion test
cch.e1.missdist <- cch.e1.miss * cch.o.e1.d

# Emotion Understanding, Threats-First Condition
cch.o.e2.d <- c()
for (i in 1:length(cch.o.e2.pca[,1])){ 
  # pull row for sub
  o.guess <- cch.o.e2.pca[i,1:pc.use]
  # transform into 1-row matrix to preserve cols
  o.guess <- matrix(o.guess,ncol = pc.use, nrow=1)
  # generate comparison hv
  o.exb <- expectation_ball(o.guess, num.samples = 5)
  # get min euclidean distance
  #set.seed(1234)
  cch.o.e2.d[i] <- hypervolume_distance(hv1 = cch.e1.hv,
                                        hv2 = o.exb,
                                        type = "minimum",
                                        num.points.max = 1000,
                                        check.memory = F)
}  

# Miss Distances
cch.e2.miss <- !cch.o.e2.inc #inverts the inclusion test
cch.e2.missdist <- cch.e2.miss * cch.o.e2.d


## C. Performance Summaries ####
imm.o.t1.res <- data.frame(ResponseID = imm.t1.o$ResponseID,
                           Danger = imm.t1.o$Danger,
                           Raw_IncTest = as.numeric(imm.o.t1.inc),
                           Raw_Dist = imm.o.t1.d, 
                           Acc_Dist = imm.t1.missdist)
head(imm.o.t1.res)
str(imm.o.t1.res)
# 108 5

cch.o.t1.res <- data.frame(ResponseID = cch.t1.o$ResponseID,
                           Danger = cch.t1.o$Danger,
                           Raw_IncTest = as.numeric(cch.o.t1.inc),
                           Raw_Dist = cch.o.t1.d,
                           Acc_Dist = cch.t1.missdist)
head(cch.o.t1.res)

imm.o.t2.res <- data.frame(ResponseID = imm.t2.o$ResponseID,
                           Danger = imm.t2.o$Danger,
                           Raw_IncTest = as.numeric(imm.o.t2.inc),
                           Raw_Dist = imm.o.t2.d,
                           Acc_Dist = imm.t2.missdist)
head(imm.o.t2.res)

cch.o.t2.res <- data.frame(ResponseID = cch.t2.o$ResponseID,
                           Danger = cch.t2.o$Danger,
                           Raw_IncTest = as.numeric(cch.o.t2.inc),
                           Raw_Dist = cch.o.t2.d,
                           Acc_Dist = cch.t2.missdist)
head(cch.o.t2.res)

imm.o.e1.res <- data.frame(ResponseID = imm.e1.o$ResponseID,
                           Danger = imm.e1.o$Danger,
                           Raw_IncTest = as.numeric(imm.o.e1.inc),
                           Raw_Dist = imm.o.e1.d,
                           Acc_Dist = imm.e1.missdist)
head(imm.o.e1.res)

cch.o.e1.res <- data.frame(ResponseID = cch.e1.o$ResponseID,
                           Danger = cch.e1.o$Danger,
                           Raw_IncTest = as.numeric(cch.o.e1.inc),
                           Raw_Dist = cch.o.e1.d,
                           Acc_Dist = cch.e1.missdist)
head(cch.o.e1.res)

imm.o.e2.res <- data.frame(ResponseID = imm.e2.o$ResponseID,
                           Danger = imm.e2.o$Danger,
                           Raw_IncTest = as.numeric(imm.o.e2.inc),
                           Raw_Dist = imm.o.e2.d,
                           Acc_Dist = imm.e2.missdist)
head(imm.o.e2.res)

cch.o.e2.res <- data.frame(ResponseID = cch.e2.o$ResponseID,
                           Danger = cch.e2.o$Danger,
                           Raw_IncTest = as.numeric(cch.o.e2.inc),
                           Raw_Dist = cch.o.e2.d,
                           Acc_Dist = cch.e2.missdist)
head(cch.o.e2.res)


## D. Binary Accuracy Correlations ####
### Threats-First Conditions ####
# Illegal Immigration
imm.o.t1e2.thr.acc <- select(imm.o.t1.res, ResponseID, Danger, Raw_IncTest) 
names(imm.o.t1e2.thr.acc) <- c("ResponseID", "Danger", "Thr_Acc")
imm.o.t1e2.emo.acc <- select(imm.o.e2.res, ResponseID, Raw_IncTest) 
names(imm.o.t1e2.emo.acc) <- c("ResponseID", "Emo_Acc")
imm.o.t1e2.all.acc <- left_join(imm.o.t1e2.thr.acc, 
                                imm.o.t1e2.emo.acc)
imm.o.t1e2.all.acc$Issue <- c("Illegal Immigration")
head(imm.o.t1e2.all.acc)

#Get correlation coefficients
imm.o.t1.all.cor <- cor(imm.o.t1e2.all.acc$Thr_Acc,
                        imm.o.t1e2.all.acc$Emo_Acc)
imm.o.t1.all.cor
#0.03432703

# Similar
imm.o.t1.sim.cor <-cor(filter(imm.o.t1e2.all.acc, Danger >= 50)$Thr_Acc,
                       filter(imm.o.t1e2.all.acc, Danger >= 50)$Emo_Acc)
imm.o.t1.sim.cor 
#-0.01138313

# Dissimilar
imm.o.t1.dis.cor <-cor(filter(imm.o.t1e2.all.acc, Danger < 50)$Thr_Acc,
                       filter(imm.o.t1e2.all.acc, Danger < 50)$Emo_Acc)
imm.o.t1.dis.cor
#0.06451903 


# Climate Change
cch.o.t1e2.thr.acc <- select(cch.o.t1.res, ResponseID, Danger, Raw_IncTest) 
names(cch.o.t1e2.thr.acc) <- c("ResponseID", "Danger","Thr_Acc")
cch.o.t1e2.emo.acc <- select(cch.o.e2.res, ResponseID, Raw_IncTest) 
names(cch.o.t1e2.emo.acc) <- c("ResponseID", "Emo_Acc")
cch.o.t1e2.all.acc <- left_join(cch.o.t1e2.thr.acc, 
                                cch.o.t1e2.emo.acc)
cch.o.t1e2.all.acc$Issue <- c("Climate Change")
head(cch.o.t1e2.all.acc)

#Get correlation coefficients
cch.o.t1.all.cor <- cor(cch.o.t1e2.all.acc$Thr_Acc,
                        cch.o.t1e2.all.acc$Emo_Acc)
cch.o.t1.all.cor
#0.2154718

cch.o.t1.sim.cor <-cor(filter(cch.o.t1e2.all.acc, Danger >= 50)$Thr_Acc,
                       filter(cch.o.t1e2.all.acc, Danger >= 50)$Emo_Acc)
cch.o.t1.sim.cor 
#0.2404588

cch.o.t1.dis.cor <-cor(filter(cch.o.t1e2.all.acc, Danger < 50)$Thr_Acc,
                       filter(cch.o.t1e2.all.acc, Danger < 50)$Emo_Acc)
cch.o.t1.dis.cor 
#0.3411211 

### Emotions-First Conditions ####
# Illegal Immigration
imm.o.t2e1.thr.acc <- select(imm.o.t2.res, ResponseID, Danger, Raw_IncTest) 
names(imm.o.t2e1.thr.acc) <- c("ResponseID", "Danger", "Thr_Acc")
imm.o.t2e1.emo.acc <- select(imm.o.e1.res, ResponseID, Raw_IncTest) 
names(imm.o.t2e1.emo.acc) <- c("ResponseID", "Emo_Acc")
imm.o.t2e1.all.acc <- left_join(imm.o.t2e1.thr.acc, 
                                imm.o.t2e1.emo.acc)
imm.o.t2e1.all.acc$Issue <- c("Illegal Immigration")
head(imm.o.t2e1.all.acc)

# Get correlation coefficients
imm.o.t2.all.cor <- cor(imm.o.t2e1.all.acc$Thr_Acc,
                        imm.o.t2e1.all.acc$Emo_Acc)
imm.o.t2.all.cor
# 0.1936701

# Similar 
imm.o.t2.sim.cor <-cor(filter(imm.o.t2e1.all.acc, Danger >= 50)$Thr_Acc,
                       filter(imm.o.t2e1.all.acc, Danger >= 50)$Emo_Acc)
imm.o.t2.sim.cor
# 0.2223474

# Dissimilar
imm.o.t2.dis.cor <-cor(filter(imm.o.t2e1.all.acc, Danger < 50)$Thr_Acc,
                       filter(imm.o.t2e1.all.acc, Danger < 50)$Emo_Acc)
imm.o.t2.dis.cor
# 0.1406937

# Climate Change
cch.o.t2e1.thr.acc <- select(cch.o.t2.res, ResponseID, Danger, Raw_IncTest) 
names(cch.o.t2e1.thr.acc) <- c("ResponseID", "Danger","Thr_Acc")
cch.o.t2e1.emo.acc <- select(cch.o.e1.res, ResponseID, Raw_IncTest) 
names(cch.o.t2e1.emo.acc) <- c("ResponseID", "Emo_Acc")
cch.o.t2e1.all.acc <- left_join(cch.o.t2e1.thr.acc, 
                                cch.o.t2e1.emo.acc)
cch.o.t2e1.all.acc$Issue <- c("Climate Change")
head(cch.o.t2e1.all.acc)

# Get correlation coefficients
cch.o.t2.all.cor <- cor(cch.o.t2e1.all.acc$Thr_Acc,
                        cch.o.t2e1.all.acc$Emo_Acc)
cch.o.t2.all.cor
# 0.2432647

# Similar 
cch.o.t2.sim.cor <-cor(filter(cch.o.t2e1.all.acc, Danger >= 50)$Thr_Acc,
                       filter(cch.o.t2e1.all.acc, Danger >= 50)$Emo_Acc)
cch.o.t2.sim.cor
# 0.3169335

# Dissimilar
cch.o.t2.dis.cor <-cor(filter(cch.o.t2e1.all.acc, Danger < 50)$Thr_Acc,
                       filter(cch.o.t2e1.all.acc, Danger < 50)$Emo_Acc)
cch.o.t2.dis.cor
# 0.04303315

#################################################

#################################################
# 5. Between-Condition Statistical Tests ####
## A. X^2 Tests of Binary Accuracy ####
### Threat Estimation Accuracy by Issue, Threats-First ####
# Immigration vs. Climate Change 
imm.o.t1.all.r1 <- table(imm.o.t1.res$Raw_IncTest)
cch.o.t1.all.r2 <- table(cch.o.t1.res$Raw_IncTest)
issue.o.t1.all.tab <- rbind(imm.o.t1.all.r1, cch.o.t1.all.r2)
colnames(issue.o.t1.all.tab) <- c("Miss", "Hit")
rownames(issue.o.t1.all.tab) <- c("Immigration", "Climate Change")

set.seed(1234)
chisq.test(issue.o.t1.all.tab, simulate.p.value = T)
# X-squared = 2.0252, df = NA, p-value = 0.1979

### Threat Estimation Accuracy by Similarity, Threats-First ####
# Threat-First: Similar vs. Dissimilar
# Immigration:
imm.o.t1.r1 <- table(filter(imm.o.t1.res, Danger >=50)$Raw_IncTest)
imm.o.t1.r2 <- table(filter(imm.o.t1.res, Danger <50)$Raw_IncTest)
imm.o.t1.tab <- rbind(imm.o.t1.r1, imm.o.t1.r2)
colnames(imm.o.t1.tab) <- c("Miss", "Hit")
rownames(imm.o.t1.tab) <- c("Similar", "Dissimilar")

set.seed(1234)
chisq.test(imm.o.t1.tab, simulate.p.value = T)
#X-squared = 5.587, df = NA, p-value = 0.02299

# Climate Change:
cch.o.t1.r1 <- table(filter(cch.o.t1.res, Danger >=50)$Raw_IncTest)
cch.o.t1.r2 <- table(filter(cch.o.t1.res, Danger <50)$Raw_IncTest)
cch.o.t1.tab <- rbind(cch.o.t1.r1, cch.o.t1.r2)
colnames(cch.o.t1.tab) <- c("Miss", "Hit")
rownames(cch.o.t1.tab) <- c("Similar", "Dissimilar")

set.seed(1234)
chisq.test(cch.o.t1.tab, simulate.p.value = T)
#X-squared = 9.2093, df = NA, p-value = 0.004498



### Threat Estimation Accuracy by Block Order, Similar Mentalizers ####
# Similar Mentalizers: Threat-First vs. Emotion-First 
imm.o.t1.r1 <- table(filter(imm.o.t1.res, Danger >=50)$Raw_IncTest)
imm.o.t2.r2 <- table(filter(imm.o.t2.res, Danger >=50)$Raw_IncTest)
imm.o.t12.tab <- rbind(imm.o.t1.r1, imm.o.t2.r2)
colnames(imm.o.t12.tab) <- c("Miss", "Hit")
rownames(imm.o.t12.tab) <- c("Threat-First", "Emotion-First")

set.seed(1234)
chisq.test(imm.o.t12.tab, simulate.p.value = T)
#X-squared = 0, df = NA, p-value = 1
# Note: these two subsets have the same N (61), which is not an error

cch.o.t1.r1 <- table(filter(cch.o.t1.res, Danger >=50)$Raw_IncTest)
cch.o.t2.r2 <- table(filter(cch.o.t2.res, Danger >=50)$Raw_IncTest)
cch.o.t12.tab <- rbind(cch.o.t1.r1, cch.o.t2.r2)
colnames(cch.o.t12.tab) <- c("Miss", "Hit")
rownames(cch.o.t12.tab) <- c("Threat-First", "Emotion-First")

set.seed(1234)
chisq.test(cch.o.t12.tab, simulate.p.value = T)
#X-squared = 0.27086, df = NA, p-value = 0.6612

### Threat Estimation Accuracy by Block Order, Dissimilar Mentalizers ####
# Dissimilar Mentalizers: Threat-First vs. Emotion-First 
imm.o.t1.r1 <- table(filter(imm.o.t1.res, Danger <50)$Raw_IncTest)
imm.o.t2.r2 <- table(filter(imm.o.t2.res, Danger <50)$Raw_IncTest)
imm.o.t12.tab <- rbind(imm.o.t1.r1, imm.o.t2.r2)
colnames(imm.o.t12.tab) <- c("Miss", "Hit")
rownames(imm.o.t12.tab) <- c("Threat-First", "Emotion-First")

set.seed(1234)
chisq.test(imm.o.t12.tab, simulate.p.value = T)
#X-squared = 2.0113, df = NA, p-value = 0.1819

cch.o.t1.r1 <- table(filter(cch.o.t1.res, Danger <50)$Raw_IncTest)
cch.o.t2.r2 <- table(filter(cch.o.t2.res, Danger <50)$Raw_IncTest)
cch.o.t12.tab <- rbind(cch.o.t1.r1, cch.o.t2.r2)
colnames(cch.o.t12.tab) <- c("Miss", "Hit")
rownames(cch.o.t12.tab) <- c("Threat-First", "Emotion-First")

set.seed(1234)
chisq.test(cch.o.t12.tab, simulate.p.value = T)
#X-squared = 0.7648, df = NA, p-value = 0.5462

## B. Mann-Whitney U Tests of Miss Distance ####
### Threat Estimation Misses by Block Order, Similar Mentalizers ####
wilcox.test(x = filter(imm.o.t1.res, Danger >=50 & Raw_IncTest == 0)$Acc_Dist,
            y = filter(imm.o.t2.res, Danger >=50 & Raw_IncTest == 0)$Acc_Dist,
            paired = F,
            exact = T)

wilcox.test(x = filter(cch.o.t1.res, Danger >=50 & Raw_IncTest == 0)$Acc_Dist,
            y = filter(cch.o.t2.res, Danger >=50 & Raw_IncTest == 0)$Acc_Dist,
            paired = F,
            exact = T)
#W = 52, p-value = 0.1832

### Threat Estimation Miss Distance by Block Order, Dissimilar Mentalizers ####
wilcox.test(x = filter(imm.o.t1.res, Danger <50 & Raw_IncTest == 0)$Acc_Dist,
            y = filter(imm.o.t2.res, Danger <50 & Raw_IncTest == 0)$Acc_Dist,
            paired = F,
            exact = T)
#W = 114, p-value = 0.7031

wilcox.test(x = filter(cch.o.t1.res, Danger <50 & Raw_IncTest == 0)$Acc_Dist,
            y = filter(cch.o.t2.res, Danger <50 & Raw_IncTest == 0)$Acc_Dist,
            paired = F,
            exact = T)
#W = 28, p-value = 0.6605
#New: W = 30, p-value = 0.8075

## C. Binary Accuracy Correlations by Block Order, Similar Mentalizers ####

# Get sample sizes
n_sim_t1_imm <-  length(filter(imm.o.t1e2.all.acc, Danger >= 50)$Thr_Acc)
n_sim_t2_imm <-  length(filter(imm.o.t2e1.all.acc, Danger >= 50)$Thr_Acc)

n_sim_t1_cch <-  length(filter(cch.o.t1e2.all.acc, Danger >= 50)$Thr_Acc)
n_sim_t2_cch <-  length(filter(cch.o.t2e1.all.acc, Danger >= 50)$Thr_Acc)

# Test whether correlations differ between block orders 
# Immigration
imm.o.sim.compcor <- cocor.indep.groups(r1.jk = imm.o.t1.sim.cor,
                                        r2.hm = imm.o.t2.sim.cor,
                                        n1 = n_sim_t1_imm,
                                        n2 = n_sim_t2_imm,
                                        alternative = "two.sided")
imm.o.sim.compcor
# fisher1925: Fisher's z (1925)
#   z = 1.2790, p-value = 0.2009
#   Null hypothesis retained
# zou2007: Zou's (2007) confidence interval
#   95% confidence interval for r1.jk - r2.hm: -0.1241 0.5720
#   Null hypothesis retained (Interval includes 0)

# Climate Change
cch.o.sim.compcor <- cocor.indep.groups(r1.jk = cch.o.t1.sim.cor,
                                        r2.hm = cch.o.t2.sim.cor,
                                        n1 = n_sim_t1_cch,
                                        n2 = n_sim_t2_cch,
                                        alternative = "two.sided")
cch.o.sim.compcor
# fisher1925: Fisher's z (1925)
#   z = -0.5246, p-value = 0.5999
#   Null hypothesis retained
# zou2007: Zou's (2007) confidence interval
#   95% confidence interval for r1.jk - r2.hm: -0.3582 0.2072
#   Null hypothesis retained (Interval includes 0)



#################################################

#################################################
# 7. Logistic Regressions on Binary Accuracy ####

## A. Prepare Data #### 
d.covars <- rbind(d.imm.oth, d.cch.oth) %>%
  select(ResponseID, Gender, Age, SelfID, PID_7, Ideo_7,
         Type, Item, PT_order, Danger)
# Recode for 3-level ID where 1 = Dem; 2 = Rep, and 3 = true Ind
d.covars$PID_3 <- as.factor(
  ifelse(as.numeric(d.covars$PID_7) < 4, 1, 
         ifelse(as.numeric(d.covars$PID_7) > 4, 2, 3)))
#table(d.covars$PID_3)/425
# Recode for 3-level ideology where 1 = Lib, 2 = Cons and 3 = Mod
d.covars$Ideo_3 <- as.factor(
  ifelse(as.numeric(d.covars$Ideo_7) < 4, 1, 
         ifelse(as.numeric(d.covars$Ideo_7) > 4, 2, 3)))
#table(d.covars$Ideo_3)/425
# Make binaries for Dem, Rep, Lib and Con
d.covars$Dem <- ifelse(d.covars$PID_3 == 1, 1, 0)
d.covars$Rep <- ifelse(d.covars$PID_3 == 2, 1, 0)
d.covars$lib <- ifelse(d.covars$Ideo_3 == 1, 1, 0)
d.covars$con <- ifelse(d.covars$Ideo_3 == 2, 1, 0)
d.covars$White <- ifelse(d.covars$SelfID == "White", 1, 0)
d.covars$ResponseID <- as.factor(d.covars$ResponseID)

# Add Task Performance
d.dists <- rbind(imm.o.t1.res, imm.o.t2.res,
                 cch.o.t1.res, cch.o.t2.res)
# IV of interest: Binary identifier for Similar Raters
d.dists$Sim_bin <- ifelse(d.dists$Danger >= 50, 1, 0)
d.mod <- left_join(d.covars, d.dists, by = c("ResponseID", "Danger"))

## B. Immigration Models ####
d.imm <- filter(d.mod, Item == "Illegalimm")

# Models Listed in Table 1:
# (1) mod.imm.sim0
# (2) mod.imm.simC
# (3) mod.imm.pty1
# (4) mod.imm.pty1C
# (5) mod.imm.ideo1
# (6) mod.imm.ideo1C
# (7) mod.imm.omni1
# (8) mod.imm.omni1C

mod.imm.sim0 <- glm(Raw_IncTest ~ Sim_bin + PT_order, 
                    data = d.imm, 
                    family = binomial(link = "logit"))
summary(mod.imm.sim0)
# bottom of the range for substantive effect:
# exp(0.7251) #2.064938 

mod.imm.simC <- update(mod.imm.sim0, . ~ . + Gender + Age + White)

mod.imm.pty1 <- glm(Raw_IncTest ~ Rep + PT_order, 
                    data = d.imm, 
                    family = binomial(link = "logit"))
mod.imm.pty1C <- update(mod.imm.pty1, . ~ . + Gender + Age + White)

mod.imm.ideo1 <- glm(Raw_IncTest ~ con + PT_order, 
                     data = d.imm, 
                     family = binomial(link = "logit"))
mod.imm.ideo1C <- update(mod.imm.ideo1, . ~ . + Gender + Age + White)

mod.imm.omni1 <- glm(Raw_IncTest ~ Sim_bin + Rep + con + 
                       PT_order, 
                     data = d.imm, 
                     family = binomial(link = "logit"))
mod.imm.omni1C <- update(mod.imm.omni1, . ~ . + Gender + Age + White)
summary(mod.imm.omni1C) 


## C. Climate Change Models ####
d.cch <- filter(d.mod, Item == "Climate")

# Models Listed in Table 1:
# (1) mod.cch.sim0
# (2) mod.cch.simC
# (3) mod.cch.pty1
# (4) mod.cch.pty1C
# (5) mod.cch.ideo1
# (6) mod.cch.ideo1C
# (7) mod.cch.omni1
# (8) mod.cch.omni1C

mod.cch.sim0 <- glm(Raw_IncTest ~ Sim_bin + PT_order, 
                    data = d.cch, 
                    family = binomial(link = "logit"))
summary(mod.cch.sim0)
# top of the range for a substantive effect
# exp(1.12693) #3.086167 - upper bound

mod.cch.simC <- update(mod.cch.sim0, . ~ . + Gender + Age + White)

mod.cch.pty1 <- glm(Raw_IncTest ~ Dem + PT_order, 
                    data = d.cch, 
                    family = binomial(link = "logit"))
mod.cch.pty1C <- update(mod.cch.pty1, . ~ . + Gender + Age + White)

mod.cch.ideo1 <- glm(Raw_IncTest ~ lib + PT_order, 
                     data = d.cch, 
                     family = binomial(link = "logit"))
mod.cch.ideo1C <- update(mod.cch.ideo1, . ~ . + Gender + Age + White)

mod.cch.omni1 <- glm(Raw_IncTest ~ Sim_bin + Dem + lib + 
                       PT_order, 
                     data = d.cch, 
                     family = binomial(link = "logit"))
mod.cch.omni1C <- update(mod.cch.omni1, . ~ . + Gender + Age + White)
summary(mod.cch.omni1C) 

#################################################

#################################################
# 8. Tests vs. Random Response Benchmarks ####

## A. Load Random Response Data (500 Datasets x 4 HVs) ####
load("AccDist_500Random_Draws_240719.RDS")

# Parse the list:
imm.t.rand.dist <- rand.hv.dist[[1]]
cch.t.rand.dist <- rand.hv.dist[[2]]
imm.e.rand.dist <- rand.hv.dist[[3]]
cch.e.rand.dist <- rand.hv.dist[[4]]
rm(rand.hv.dist)

## B. Derive Random Response Benchmarks ####
### Binary Accuracy: Threat Estimation ####

# Get the binary accuracy rate for each of the 500 samples
imm.t.rand.acc <- sapply(imm.t.rand.dist, function(x) mean(x[,2]))
cch.t.rand.acc <- sapply(cch.t.rand.dist, function(x) mean(x[,2]))

# Make dataframes for plotting later
imm.t.rand.acc.df <- data.frame(Iteration = c(1:500),
                                 Accuracy = imm.t.rand.acc)
cch.t.rand.acc.df <- data.frame(Iteration = c(1:500),
                                 Accuracy = cch.t.rand.acc)

### Binary Accuracy: Emotion Understanding ####
imm.e.rand.acc <- sapply(imm.e.rand.dist, function(x) mean(x[,2]))
cch.e.rand.acc <- sapply(cch.e.rand.dist, function(x) mean(x[,2]))

imm.e.rand.acc.df <- data.frame(Iteration = c(1:500),
                                 Accuracy = imm.e.rand.acc)
cch.e.rand.acc.df <- data.frame(Iteration = c(1:500),
                                 Accuracy = cch.e.rand.acc)

### Miss Distance: Threat Estimation ####
imm.t.rand.dmiss <- lapply(imm.t.rand.dist, function(x) filter(x, Guess_dist > 0))
cch.t.rand.dmiss <- lapply(cch.t.rand.dist, function(x) filter(x, Guess_dist > 0))

# Check distributions to determine appropriate summary stat
shapiro.test(imm.t.rand.dmiss[[1]]$Guess_dist)
#W = 0.8903, p-value = 0.0003582
shapiro.test(cch.t.rand.dmiss[[1]]$Guess_dist)
#W = 0.89492, p-value = 0.001018 
# Assume non-normal; use medians

# Get median miss distances
imm.t.rand.dmiss.med <- sapply(imm.t.rand.dmiss, function(x) median(x[,"Guess_dist"]))
cch.t.rand.dmiss.med <- sapply(cch.t.rand.dmiss, function(x) median(x[,"Guess_dist"]))

# Make dataframes for plotting
imm.t.rand.dmiss.df <- data.frame(Iteration = c(1:500),
                                   Median_Miss = imm.t.rand.dmiss.med)
cch.t.rand.dmiss.df <- data.frame(Iteration = c(1:500),
                                   Median_Miss = cch.t.rand.dmiss.med)

### Correlation of Binary Accuracies ####

# Merge lists to merge dfs by iteration
imm.rand.all <- purrr::map2(imm.t.rand.dist,
                            imm.e.rand.dist,
                            inner_join, 
                            by = "Obs")
cch.rand.all <- purrr::map2(cch.t.rand.dist,
                            cch.e.rand.dist,
                            inner_join, 
                            by = "Obs")

# Get Correlations
# note: .x is threats, .y is emotions
imm.rand.corr <- sapply(imm.rand.all, function(x) cor(x[,"Inc_test.x"],
                                                      x[,"Inc_test.y"]))
cch.rand.corr <- sapply(cch.rand.all, function(x) cor(x[,"Inc_test.x"],
                                                      x[,"Inc_test.y"]))

## C. True Binary Accuracy vs. Benchmark ####

# Immigration:
# Get accuracy rate from real data
imm.t1.acc.true <- mean(imm.o.t1.res$Raw_IncTest)
imm.t1.acc.true
#0.712963
# Compare accuracy rate to chance distribution (non-parametric exact p)
sum(imm.t.rand.acc >= imm.t1.acc.true) 
# 5/500: p = 0.01

# Partition by Similiarty
imm.t1.acc.sim.true <- mean(filter(imm.o.t1.res, Danger >= 50)$Raw_IncTest)
imm.t1.acc.sim.true
#0.8032787
imm.t1.acc.dis.true <- mean(filter(imm.o.t1.res, Danger < 50)$Raw_IncTest)
imm.t1.acc.dis.true
#0.5957447
# Compare True Accuracy to Random Draw (p-value of 0.05 is equiv to < 25)
sum(imm.t.rand.acc >= imm.t1.acc.sim.true) 
#0/500, p < 0.002
sum(imm.t.rand.acc >= imm.t1.acc.dis.true) 
#237/500, p = 0.474

# Climate Change:
cch.t1.acc.true <- mean(cch.o.t1.res$Raw_IncTest)
cch.t1.acc.true
#0.7962963
sum(cch.t.rand.acc >= cch.t1.acc.true) #0
# 0/500: p < 0.002 (sensitivity limit of the test)

cch.t1.acc.sim.true <- mean(filter(cch.o.t1.res, Danger >= 50)$Raw_IncTest)
cch.t1.acc.sim.true
#0.8641975
cch.t1.acc.dis.true <- mean(filter(cch.o.t1.res, Danger < 50)$Raw_IncTest)
cch.t1.acc.dis.true
#0.5925926
sum(cch.t.rand.acc >= cch.t1.acc.sim.true) 
#0/500, p < 0.002
sum(cch.t.rand.acc >= cch.t1.acc.dis.true) 
#105/500, p = 0.21

## D. True Miss Distance vs. Benchmark ####
# Get the true misses 
imm.t1.dmiss.true <- filter(imm.o.t1.res, Acc_Dist > 0)
cch.t1.dmiss.true <- filter(cch.o.t1.res, Acc_Dist > 0)

# Check distribution - not normal, use medians
shapiro.test(imm.t1.dmiss.true$Acc_Dist) 
#W = 0.72415, p-value = 2.694e-06
shapiro.test(cch.t1.dmiss.true$Acc_Dist) 
#W = 0.89355, p-value = 0.02213

# Immigration Threat Estimation:
imm.t1.dmiss.true.med <- median(imm.t1.dmiss.true$Acc_Dist)
imm.t1.dmiss.true.med
# 16.67041
# Guesses are better if distances are smaller
sum(imm.t.rand.dmiss.med <= imm.t1.dmiss.true.med) 
#1/500, p = 0.002

# Partition by Similiarty
imm.t1.dmiss.sim.true.med <- median(filter(imm.t1.dmiss.true, Danger >= 50)$Acc_Dist)
imm.t1.dmiss.dis.true.med <- median(filter(imm.t1.dmiss.true, Danger < 50)$Acc_Dist)
# Compare Distance of Misses to Random Draw 
sum(imm.t.rand.dmiss.med <= imm.t1.dmiss.sim.true.med) 
#0/500, p < 0.002
sum(imm.t.rand.dmiss.med <= imm.t1.dmiss.dis.true.med) 
#6/500, p = 0.012

# Climate Change Threat Estimation:
cch.t1.dmiss.true.med <- median(cch.t1.dmiss.true$Acc_Dist)
sum(cch.t.rand.dmiss.med <= cch.t1.dmiss.true.med) 
#0/500, p < 0.002

# Partition by Similiarty
cch.t1.dmiss.sim.true.med <- median(filter(cch.t1.dmiss.true, Danger >= 50)$Acc_Dist)
cch.t1.dmiss.dis.true.med <- median(filter(cch.t1.dmiss.true, Danger < 50)$Acc_Dist)
# Compare Distance of Misses to Random Draw
sum(cch.t.rand.dmiss.med <= cch.t1.dmiss.sim.true.med) 
#0/500, p < 0.002
sum(cch.t.rand.dmiss.med <= cch.t1.dmiss.dis.true.med) 
#1/500, p = 0.002

## E. True Accuracy Correlations vs. Benchmark  ####
# Emotions-First Conditions
# All Mentalizers
sum(imm.rand.corr >= imm.o.t2.all.cor) 
#28/500, p = 0.056
sum(cch.rand.corr >= cch.o.t2.all.cor) 
#34/500, p = 0.068

# Similar
sum(imm.rand.corr >= imm.o.t2.sim.cor) 
#12/500, p = 0.024
sum(cch.rand.corr >= cch.o.t2.sim.cor) 
#7/500, p = 0.014

# Dissimilar
sum(imm.rand.corr >= imm.o.t2.dis.cor) 
#73/500, p = 0.146
sum(cch.rand.corr >= cch.o.t2.dis.cor) 
#361/500, p = 0.722



#################################################

#################################################
# 9. Main Text Figures ####

# NOTE: compound figures created in Adobe Illustrator
# Hypervolume snapshots for Figure 2 saved as pngs
# Input images for Figures 3-5 saved as pdfs set to 5x7 in scale

## Figure 2A: Hypervolume for Illegal Immigration Threat Perception ####

# Redefine Mentalizer data for plotting
imm.o.t1.plot.df <- imm.o.t1.res

# Dataframe with xyz values and random vs. other points as factor
imm.t1.hv.rps <- data.frame(imm.t1.hv@RandomPoints)
imm.t1.hv.rps$PType <- c("Random")
# Get ground truth values
imm.t1.hv.grd <- data.frame(imm.t1.hv@Data)
imm.t1.hv.grd$PType <- c("Self-Rating")

# Add coords for guesses from the Other ratings
imm.o.t1.plot.df[1:20,1:5]
# Sub 1 = accurate guess
# Sub 8 = good guess
# Sub 12 = bad guess
imm.t1.hv.guesses <- rbind(imm.o.t1.pca[1,1:pc.use],
                           imm.o.t1.pca[8,1:pc.use],
                           imm.o.t1.pca[12,1:pc.use])
imm.t1.hv.guesses <- data.frame(imm.t1.hv.guesses, 
                                PType = c("Accurate Guess",
                                          "Good Guess",
                                          "Bad Guess"))
# For Plotting
imm.t1.hv.toplot <- rbind(imm.t1.hv.rps,
                          imm.t1.hv.grd,
                          imm.t1.hv.guesses)
imm.t1.hv.toplot$PType <- factor(imm.t1.hv.toplot$PType)
imm.t1.hv.toplot$Colors <- ifelse(imm.t1.hv.toplot$PType == "Random", "#cccccc",
                                  ifelse(imm.t1.hv.toplot$PType == "Self-Rating", "#252525",
                                         "#de2d26"))
imm.t1.hv.toplot$Colors <- factor(imm.t1.hv.toplot$Colors)

# Plot Settings
dev.off()
plot3d(x = imm.t1.hv.toplot$PC1,
       y = imm.t1.hv.toplot$PC2,
       z = imm.t1.hv.toplot$PC3,
       xlab = "PC1", ylab = "PC2", zlab = "PC3",
       type = "s", radius=2, 
       col = imm.t1.hv.toplot$Colors)
# Play with rotation
# rgl.snapshot('3dPlot_imm.png', fmt='png')

## Figure 2B: Hypervolume for Climate Change Threat Perception ####

# Redefine Mentalizer data for plotting
cch.o.t1.plot.df <- cch.o.t1.res

# Create df. with xyz values and random vs. other points as factor
cch.t1.hv.rps <- data.frame(cch.t1.hv@RandomPoints)
cch.t1.hv.rps$PType <- c("Random")

# Get ground truth values
cch.t1.hv.grd <- data.frame(cch.t1.hv@Data)
cch.t1.hv.grd$PType <- c("Self-Rating")

# Add coords for guesses from the Other ratings
cch.o.t1.plot.df[1:20,1:5]
# Sub 1 = accurate guess
# Sub 8 = good guess
# Sub 15 = bad guess

cch.t1.hv.guesses <- rbind(cch.o.t1.pca[1,1:pc.use],
                           cch.o.t1.pca[8,1:pc.use],
                           cch.o.t1.pca[15,1:pc.use])
cch.t1.hv.guesses <- data.frame(cch.t1.hv.guesses, 
                                PType = c("Accurate Guess",
                                          "Good Guess",
                                          "Bad Guess"))
# For Plotting
cch.t1.hv.toplot <- rbind(cch.t1.hv.rps,
                          cch.t1.hv.grd,
                          cch.t1.hv.guesses)
cch.t1.hv.toplot$PType <- factor(cch.t1.hv.toplot$PType)
cch.t1.hv.toplot$Colors <- ifelse(cch.t1.hv.toplot$PType == "Random", "#cccccc",
                                  ifelse(cch.t1.hv.toplot$PType == "Self-Rating", "#252525",
                                         "#de2d26"))
cch.t1.hv.toplot$Colors <- factor(cch.t1.hv.toplot$Colors)

# Plot 
dev.off()
plot3d(x = cch.t1.hv.toplot$PC1,
       y = cch.t1.hv.toplot$PC2,
       z = cch.t1.hv.toplot$PC3,
       xlab = "PC1", ylab = "PC2", zlab = "PC3",
       type = "s", radius=2, 
       col = cch.t1.hv.toplot$Colors)
#rgl.snapshot('3dPlot_cch.png', fmt='png')


## Figure 3A: Density Plot for II Binary vs Chance ####

# vlines to show
v.imm.all.tru <- imm.t1.acc.true
v.imm.sim.tru <- imm.t1.acc.sim.true  
v.imm.dis.tru <- imm.t1.acc.dis.true

# Create plot w/ vlines for observed values
fig3.imm <- ggplot(imm.t.rand.acc.df,
                   aes(x = Accuracy*100)) +
  geom_histogram(fill = "#969696", 
              alpha = 0.8,
              binwidth = 2) +
  geom_vline(xintercept = v.imm.all.tru*100, 
             size = 1.5) +
  geom_vline(xintercept = v.imm.sim.tru*100, 
             size = 1.5,
             linetype = 2) +
  geom_vline(xintercept = v.imm.dis.tru*100, 
             size = 1.5,
             linetype = 3) +
  annotate("text", 
           x = 70.5, y = 84,
           label = "All\n Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 79.5, y = 84,
           label = "Similar\n Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 59, y = 84,
           label = "Dissimilar\n Mentalizers",
           hjust = 1,
           size = 5) + 
  labs(title = "Illegal Immigration: Binary Accuracy",
       subtitle = "Chance distribution based on 500 simulations",
       x = "Accuracy (% Guesses inside Ground Truth)",
       y = "Count") +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 16))
fig3.imm




## Figure 3B: Density Plot for CC Binary vs Chance ####

# vlines to show
v.cch.all.tru <- cch.t1.acc.true
v.cch.sim.tru <- cch.t1.acc.sim.true  
v.cch.dis.tru <- cch.t1.acc.dis.true

# Create plot w/ vlines for observed values
fig3.cch <- ggplot(cch.t.rand.acc.df,
                   aes(x = Accuracy*100)) +
  geom_histogram(fill = "#969696", 
                 alpha = 0.8,
                 binwidth = 2) +
  geom_vline(xintercept = v.imm.all.tru*100, 
             size = 1.5) +
  geom_vline(xintercept = v.imm.sim.tru*100, 
             size = 1.5,
             linetype = 2) +
  geom_vline(xintercept = v.imm.dis.tru*100, 
             size = 1.5,
             linetype = 3) +
  annotate("text", 
           x = 70.5, y = 80,
           label = "All\n Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 79.5, y = 80,
           label = "Similar\n Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 59, y = 80,
           label = "Dissimilar\n Mentalizers",
           hjust = 1,
           size = 5) + 
  labs(title = "Climate Change: Binary Accuracy",
       subtitle = "Chance distribution based on 500 simulations",
       x = "Accuracy (% Guesses inside Ground Truth)",
       y = "Count") +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 16))
fig3.cch



## Figure 3C: Density Plot for II: Miss Distance vs Chance ####

# vlines to show
v.imm.all.med.tru <- imm.t1.dmiss.true.med
v.imm.sim.med.tru <- imm.t1.dmiss.sim.true.med
v.imm.dis.med.tru <- imm.t1.dmiss.dis.true.med

# Create plot w/ vlines for observed values
fig3c.imm <- ggplot(imm.t.rand.dmiss.df,
                   aes(x = Median_Miss)) +
  geom_histogram(fill = "#969696", 
                 alpha = 0.8,
                 binwidth = 1) +
  geom_vline(xintercept = v.imm.all.med.tru, 
             size = 1.5) +
  geom_vline(xintercept = v.imm.sim.med.tru, 
             size = 1.5,
             linetype = 2) +
  geom_vline(xintercept = v.imm.dis.med.tru, 
             size = 1.5,
             linetype = 3) +
  annotate("text", 
           x = 16.4, y = 110,
           label = "All Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 15, y = 90,
           label = "Similar\n Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 18.5, y = 90,
           label = "Dissimilar\nMentalizers",
           hjust = 0,
           size = 5) + 
  labs(title = "Illegal Immigration: Miss Distance",
       subtitle = "Chance distribution based on 500 simulations",
       x = "Median Distance to Ground Truth",
       y = "Count") +
  xlim(10,35) +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 16))
fig3c.imm

## Figure 3D: Density Plot for CC: Miss Distance vs Chance ####

# vlines to show
v.cch.all.med.tru <- cch.t1.dmiss.true.med
v.cch.sim.med.tru <- cch.t1.dmiss.sim.true.med
v.cch.dis.med.tru <- cch.t1.dmiss.dis.true.med

# Create plot w/ vlines for observed values
fig3d.cch <- ggplot(cch.t.rand.dmiss.df,
                    aes(x = Median_Miss)) +
  geom_histogram(fill = "#969696", 
                 alpha = 0.8,
                 binwidth = 1) +
  geom_vline(xintercept = v.cch.all.med.tru, 
             size = 1.5) +
  geom_vline(xintercept = v.cch.sim.med.tru, 
             size = 1.5,
             linetype = 2) +
  geom_vline(xintercept = v.cch.dis.med.tru, 
             size = 1.5,
             linetype = 3) +
  annotate("text", 
           x = 15.3, y = 95,
           label = "All\n Mentalizers",
           hjust = 1,
           size = 5) + 
  annotate("text", 
           x = 10.8, y = 70,
           label = "Similar\nMentalizers",
           hjust = 0,
           size = 5) + 
  annotate("text", 
           x = 19.5, y = 70,
           label = "Dissimilar\nMentalizers",
           hjust = 0,
           size = 5) + 
  labs(title = "Climate Change: Miss Distance",
       subtitle = "Chance distribution based on 500 simulations",
       x = "Median Distance to Ground Truth",
       y = "Count") +
  ylim(0,100) +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 16))
fig3d.cch


## Figure 4A: Bar Plot for II by Condition  ####

imm.o.cond.st <- table(filter(imm.o.t1.res, Danger >=50)$Raw_IncTest)
imm.o.cond.se <- table(filter(imm.o.t2.res, Danger >=50)$Raw_IncTest)
imm.o.cond.dt <- table(filter(imm.o.t1.res, Danger <50)$Raw_IncTest)
imm.o.cond.de <- table(filter(imm.o.t2.res, Danger <50)$Raw_IncTest)

imm.o.cond.df <- data.frame(Condition = c(rep("Threats-First", 2),
                                          rep("Emotions-First", 2)),
                            Similarity = c("Similar Mentalizers",
                                           "Dissimilar Mentalizers",
                                           "Similar Mentalizers",
                                           "Dissimilar Mentalizers"),
                            Share = c(imm.o.cond.st[2]/sum(imm.o.cond.st),
                                      imm.o.cond.dt[2]/sum(imm.o.cond.dt),
                                      imm.o.cond.se[2]/sum(imm.o.cond.se),
                                      imm.o.cond.de[2]/sum(imm.o.cond.de)))

imm.o.cond.df$Condition <- factor(imm.o.cond.df$Condition,
                                  levels = c("Threats-First",
                                             "Emotions-First"),
                                  labels = c("Threats-\nFirst",
                                             "Emotions-\nFirst"))
imm.o.cond.df$Similarity <- factor(imm.o.cond.df$Similarity,
                                  levels = c("Similar Mentalizers",
                                             "Dissimilar Mentalizers"))

# Chance with probs 95: 0.6759259
# But see notes above on distribution
imm.chance <- quantile(imm.t.rand.acc, probs = 0.95)

fig4a.imm <- ggplot(imm.o.cond.df,
                    aes(x = Condition,
                        y = Share*100,
                        fill = Condition)) +
  facet_grid(. ~ Similarity) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = imm.chance * 100,
             linetype = "dashed") +
  scale_fill_grey() +
  labs(title = "Illegal Immigration: Binary Accuracy \nby Condition",
       y = "Accuracy (%)") +
  annotate("text", 
           x = 1.5, y = 95,
           label = "N.S.",
           size = 5) +
  geom_segment(x = 1, y = 90, xend = 2, yend = 90,
               size = 1) +
  ylim(0,100) +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.title.x = element_blank(),
        axis.text = element_text(size = 16),
        strip.text = element_text(size = 16))
fig4a.imm

## Figure 4B: Bar Plot for CC by Condition  ####

cch.o.cond.st <- table(filter(cch.o.t1.res, Danger >=50)$Raw_IncTest)
cch.o.cond.se <- table(filter(cch.o.t2.res, Danger >=50)$Raw_IncTest)
cch.o.cond.dt <- table(filter(cch.o.t1.res, Danger <50)$Raw_IncTest)
cch.o.cond.de <- table(filter(cch.o.t2.res, Danger <50)$Raw_IncTest)

cch.o.cond.df <- data.frame(Condition = c(rep("Threats-First", 2),
                                          rep("Emotions-First", 2)),
                            Similarity = c("Similar Mentalizers",
                                           "Dissimilar Mentalizers",
                                           "Similar Mentalizers",
                                           "Dissimilar Mentalizers"),
                            Share = c(cch.o.cond.st[2]/sum(cch.o.cond.st),
                                      cch.o.cond.dt[2]/sum(cch.o.cond.dt),
                                      cch.o.cond.se[2]/sum(cch.o.cond.se),
                                      cch.o.cond.de[2]/sum(cch.o.cond.de)))

cch.o.cond.df$Condition <- factor(cch.o.cond.df$Condition,
                                  levels = c("Threats-First",
                                             "Emotions-First"),
                                  labels = c("Threats-\nFirst",
                                             "Emotions-\nFirst"))
cch.o.cond.df$Similarity <- factor(cch.o.cond.df$Similarity,
                                   levels = c("Similar Mentalizers",
                                              "Dissimilar Mentalizers"))
# Chance w 5% false positives = 0.6203704
cch.chance <- quantile(cch.t.rand.acc, probs = 0.95)


fig4b.cch <- ggplot(cch.o.cond.df,
                    aes(x = Condition,
                        y = Share*100,
                        fill = Condition)) +
  facet_grid(. ~ Similarity) +
  geom_bar(stat = "identity") +
  geom_hline(yintercept = cch.chance * 100,
             linetype = "dashed") +
  scale_fill_grey() +
  labs(title = "Climate Change: Binary Accuracy \nby Condition",
       y = "Accuracy (%)") +
  annotate("text", 
           x = 1.5, y = 95,
           label = "N.S.",
           size = 5) +
  geom_segment(x = 1, y = 90, xend = 2, yend = 90,
               size = 1) +
  ylim(0,100) +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.title.x = element_blank(),
        axis.text = element_text(size = 16),
        strip.text = element_text(size = 16))
fig4b.cch


## Figure 5A: Density Plot for II Correlations ####

# vlines to show
imm.o.t2.all.cor 
imm.o.t2.sim.cor
imm.o.t2.dis.cor

imm.t2.rand.corr.df <- data.frame(Iteration = c(1:500),
                                  Correlation = imm.rand.corr)
#head(imm.t2.rand.corr.df)

# Create plot w/ vlines for observed values
fig5.imm <- ggplot(imm.t2.rand.corr.df,
                   aes(x = Correlation)) +
  geom_histogram(fill = "#969696", 
                 alpha = 0.8,
                 binwidth = 0.02) +
  geom_vline(xintercept = imm.o.t2.all.cor, 
             size = 1.5) +
  geom_vline(xintercept = imm.o.t2.sim.cor, 
             size = 1.5,
             linetype = 2) +
  geom_vline(xintercept = imm.o.t2.dis.cor, 
             size = 1.5,
             linetype = 3) +
  annotate("text",
           x = 0.185, y = 45,
           label = "All\n Mentalizers",
           hjust = 1,
           size = 5) +
  annotate("text",
           x = 0.23, y = 30,
           label = "Similar\nMentalizers",
           hjust = 0,
           size = 5) +
  annotate("text",
           x = 0.13, y = 30,
           label = "Dissimilar\nMentalizers",
           hjust = 1,
           size = 5) +
  labs(title = "Illegal Immigration: Accuracy Correlations \nin the Emotions-First Condition",
       subtitle = "Chance distribution based on 500 simulations",
       x = "Correlation (phi) between Accuracy in \nEmotion Understanding and Threat Estimation",
       y = "Count") +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 16))
fig5.imm


## Figure 5B: Density Plot for CC Correlations ####

# vlines to show
cch.o.t2.all.cor
cch.o.t2.sim.cor
cch.o.t2.dis.cor

# Create plot w/ vlines for observed values
cch.t2.rand.corr.df <- data.frame(Iteration = c(1:500),
                                  Correlation = cch.rand.corr)

# Create plot w/ vlines for observed values
fig5.cch <- ggplot(cch.t2.rand.corr.df,
                   aes(x = Correlation)) +
  geom_histogram(fill = "#969696", 
                 alpha = 0.8,
                 binwidth = 0.02) +
  geom_vline(xintercept = cch.o.t2.all.cor, 
             size = 1.5) +
  geom_vline(xintercept = cch.o.t2.sim.cor, 
             size = 1.5,
             linetype = 2) +
  geom_vline(xintercept = cch.o.t2.dis.cor, 
             size = 1.5,
             linetype = 3) +
  annotate("text",
           x = 0.23, y = 45,
           label = "All\n Mentalizers",
           hjust = 1,
           size = 5) +
  annotate("text",
           x = 0.33, y = 30,
           label = "Similar\nMentalizers",
           hjust = 0,
           size = 5) +
  annotate("text",
           x = 0.03, y = 30,
           label = "Dissimilar\nMentalizers",
           hjust = 1,
           size = 5) +
  labs(title = "Climate Change: Accuracy Correlations \nin the Emotions-First Condition",
       subtitle = "Chance distribution based on 500 simulations",
       x = "Correlation (phi) between Accuracy in \nEmotion Understanding and Threat Estimation",
       y = "Count") +
  theme_bw() +
  theme(legend.position = "none",
        title = element_text(size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 16))
fig5.cch





#################################################

#################################################
# 10. Main Text Tables ####
## Table 1: Models of Binary Threat Estimation Accuracy: II ####

# Covariate labels
imm.covars.full <- c("Similarity", "Republican", "Conservative", 
                 "Emotion-First", 
                 "Female", "Age", "White",
                 "Constant")


# toggle 'type' and specify 'out' to produce exportable table 

imm.glm.tab <- stargazer::stargazer(
  mod.imm.sim0, mod.imm.simC, mod.imm.pty1,
  mod.imm.pty1C, mod.imm.ideo1, mod.imm.ideo1C,
  mod.imm.omni1, mod.imm.omni1C,
  type="text",
#  type="html",
#  out = "imm_mods_tab.html",
  dep.var.caption = 'Issue: Illegal Immigration',
  dep.var.labels = "",
  model.numbers = T,
  no.space = T,
  digits=2,
  single.row = F,
  star.cutoffs = c(0.05, 0.01, 0.001),
  covariate.labels = imm.covars.full,
  column.sep.width = "-9pt"
)

## Table 2: Models of Binary Threat Estimation Accuracy: CC ####
cch.covars.full <- c("Similarity", "Democrat", "Liberal", 
                 "Emotion-First", 
                 "Female", "Age", "White",
                 "Constant")

# toggle 'type' and specify 'out' to produce exportable table 

stargazer::stargazer(
  mod.cch.sim0, mod.cch.simC, mod.cch.pty1,
  mod.cch.pty1C, mod.cch.ideo1, mod.cch.ideo1C,
  mod.cch.omni1, mod.cch.omni1C,
  type="text",
#  type="html",
#  out = "cch_mods_tab.html",
  dep.var.caption = 'Issue: Climate Change',
  dep.var.labels = "",
  model.numbers = T,
  no.space = T,
  digits=2,
  single.row = F,
  star.cutoffs = c(0.05, 0.01, 0.001),
  covariate.labels = cch.covars.full,
  column.sep.width = "-9pt"
)


#################################################

#################################################
# 11. Appendix Tables ####
## Table A1: Sample Demographics by Condition ####


## Stats for the Full Sample:
all.demos <- d.expt %>% 
  mutate(w = ifelse(SelfID == "White",1,0),
         college = ifelse(Edu %in% c("3", "4", "5"), 1, 0),
         dem = ifelse(PID_7 %in% c("1", "2"), 1, 0),
         indep = ifelse(PID_7 == c("4"), 1, 0),
         rep = ifelse(PID_7 %in% c("6", "7"), 1, 0),
         lib = ifelse(Ideo_7 %in% c("1", "2"), 1, 0),
         mod = ifelse(Ideo_7 %in% c("4"), 1, 0),
         cons = ifelse(Ideo_7 %in% c("6", "7"), 1, 0)) %>%
  summarize(Condition = "Full Sample",
            N = n(),
            f.share = 100* mean(as.numeric(Gender)-1),
            age.mean = mean(Age),
            age.sd = sd(Age),
            w.share = 100* mean(as.numeric(w)),
            coll.share = 100* mean(as.numeric(college)),
            dem.share = 100* mean(as.numeric(dem)),
            ind.share = 100* mean(as.numeric(indep)),
            rep.share = 100* mean(as.numeric(rep)),
            lib.share = 100* mean(as.numeric(lib)),
            mod.share = 100* mean(as.numeric(mod)),
            cons.share = 100* mean(as.numeric(cons)))
all.demos

## Stats by Condition
cond.demos <- d.expt %>% 
  mutate(w = ifelse(SelfID == "White",1,0),
         college = ifelse(Edu %in% c("3", "4", "5"), 1, 0),
         dem = ifelse(PID_7 %in% c("1", "2"), 1, 0),
         indep = ifelse(PID_7 == c("4"), 1, 0),
         rep = ifelse(PID_7 %in% c("6", "7"), 1, 0),
         lib = ifelse(Ideo_7 %in% c("1", "2"), 1, 0),
         mod = ifelse(Ideo_7 %in% c("4"), 1, 0),
         cons = ifelse(Ideo_7 %in% c("6", "7"), 1, 0)) %>%
  group_by(Condition) %>%
  summarize(N = n(),
            f.share = 100* mean(as.numeric(Gender)-1),
            age.mean = mean(Age),
            age.sd = sd(Age),
            w.share = 100* mean(as.numeric(w)),
            coll.share = 100* mean(as.numeric(college)),
            dem.share = 100* mean(as.numeric(dem)),
            ind.share = 100* mean(as.numeric(indep)),
            rep.share = 100* mean(as.numeric(rep)),
            lib.share = 100* mean(as.numeric(lib)),
            mod.share = 100* mean(as.numeric(mod)),
            cons.share = 100* mean(as.numeric(cons)))

# Combined Table:
demos.mat <- rbind(all.demos, cond.demos)
demos.df <- t(data.frame(demos.mat[,-1]))
demos.df <- round(demos.df, 1)
colnames(demos.df) <- unname(t(demos.mat[,1]))
# Produce the table
print(demos.df)

# Additional reported demographics (full sample):

# % Black or African American; % Latino or Hispanic
table(d.expt$SelfID)/length(d.expt$SelfID) 
# Asian      Black        DNR   Lat/Hisp      Mixed      Other      White 
# 0.04410012 0.12395709 0.00000000 0.10250298 0.04052443 0.02860548 0.66030989 

# % 7-point Party ID 
table(d.expt$PID_7)/length(d.expt$PID_7)
# 1          2          3          4          5          6          7 
# 0.23837902 0.15017878 0.09654350 0.14898689 0.06436234 0.11799762 0.17878427 
# 99 
# 0.00476758 
# 1 = strong democrat
# 2 = weak dem
# 3 = indep (lean dem)
# 4 = indep 
# 5 = indep (lean rep)
# 6 = weak rep
# 7 = strong republican

## Table A3: Threat Perception First Three Principal Components ####

# Identifies PC1, PC2, and PC3 loadings
# Climate Change: Columns 1-3
round(cch.t1.pca$rotation[,1:3], 2)
# Illegal Immigration: Columns 4-6 
round(imm.t1.pca$rotation[,1:3], 2)

## Table A4: Emotion Response First Three Principal Components ####

# Identifies PC1, PC2, and PC3 loadings
# Climate Change: Columns 1-3
round(cch.e1.pca$rotation[,1:3], 2)
# Illegal Immigration: Columns 4-6 
round(imm.e1.pca$rotation[,1:3], 2)



